/*
 * 单原子自旋可以向上、向下，磁矩为mu，外场B向上，哈密顿量为
 *  H = - mu B Sz
 *  下面的代码输出马尔科夫链在屏幕上
 */

#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>

void markov_single_atom(int n, double muB, double kT, double & ave_E, double & ave_Sz){

	double sum_E = 0, sum_Sz = 0;
	int Sz = 1;
	for(int i=0;i<n;i++){
		if(Sz ==1){
			//cout<<"↑";
			sum_E += - muB;
			sum_Sz += Sz;
			if( (double)rand()/RAND_MAX  < exp(-2*muB/kT) )
				Sz = -1;
		}
		else{
			sum_E += muB;
			sum_Sz += Sz;
			//cout<<"↓";
			Sz = 1;
		}
	}
	cout<<endl;

	ave_E = sum_E /n;
	ave_Sz = sum_Sz /n;
}

int main(){

	int n=10000;
	double muB = 2;
	double kT, ave_E, ave_Sz;
	ofstream fp("mc_markov_demon_sampling.txt");
	for(kT = 0; kT <50; kT += 0.1){
		markov_single_atom(n, muB, kT, ave_E, ave_Sz);
		fp<< kT<<"\t"<<ave_E<<"\t"<<ave_Sz<<endl;
	}
	fp.close();
	return 0;
}
